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1. Abstract 

Recently, a new connection between density functional theory and kinetic theory has been 
proposed. In particular, it was shown that the Kohn-Sham (KS) equations can be reformulated 
as a macroscopic limit of the steady-state solution of a suitable single-particle kinetic equation. 
By using a discrete version of this new formalism, the exchange and correlation energies of simple 
atoms and the geometrical configuration of the methane molecule were calculated accurately. 
Here, we discuss the main ideas behind the lattice kinetic approach to electronic structure 
computations, offer some considerations for prospective extensions, and also show additional 
numerical results, namely the geometrical configuration of the water molecule. 

2. Introduction 

The calculation of physical properties of interacting many-body quantum systems is one of the 
major challenges in chemistry and condensed matter. In principle, this requires the solution 
of the Schrodinger or Dirac equations for 3 N spatial variables and N spin electronic variables, 
where N is the number of particles in the system. Usually, for atoms, the number of electrons 
is in the range N ~ 1 — 100, for small molecules often more than 100, and for solids, one has 
around N ~ 10 2 ' 3 . Therefore, developing approximate models to describe these systems becomes 
crucial. In many cases, the most important goal is to calculate several measurable quantities, 
such as the bonding energy, polarisability, conductivity, etc., rather than the wave function 
itself. A formalism which can treat the systems mentioned with way less computational effort, 
while still being formally exact, is the density functional method |T], developed by Hohenberg 
and Kohn [2| and Kohn and Sham j3]. The Kohn-Sham approach to density functional theory 
allows an exact description of the interacting many-particle systems in terms of an effective non¬ 
interacting particle system. The effective potential can be shown to be completely determined 
by the total electron density of the interacting system, which is the reason why it is called a 
density functional. In particular, the ground state energy of the system is a density functional, 
whose exact expression is not known due to the complicated nature of the quantum many-body 
problem. However over the years, thanks to an intense and thoughtful work, many increasingly 
more accurate approximations have continued to appear 110 . 


On an apparently very different front, fluid dynamics, the Lattice Boltzmann (LB) method 
has been used with considerable success to simulate a variety of physical systems [6i [TJ [8]. In 
the recent years, however, it has become increasingly apparent that the LB methodology is not 
confined to fluid dynamics and can indeed be extended to a variety of complex dynamical systems 
described by linear and non-linear partial differential equations, such as electromagnetism, 
quantum mechanics, Bose-Einstein condensation and many other mmm- In a very recent 
paper, Mendoza et al. [12] showed for the first time that the LB extends to the case of electronic 
structure simulations; more precisely, these authors developed and validated a lattice kinetic 
formulation of the Kohn-Sham equations of electrons density functional theory. In this paper, 
we review the main ideas and motivations behind the above formulation and discuss main 
challenges ahead. 


3. Electron density functional theory: Kohn-Sham equations 

The quantum mechanical behaviour of a molecule consisting of Nj ions and N e electrons is 
described by the Schrodinger equation for the N-body wave function ^/(ri... rjv e ; R\ ■ ■ ■ RnR t), 
that is 

ihdpV = HV , (1) 

where H is the N-body Hamiltonian of the system, i.e. kinetic energy plus electron-electron and 
electron-ion Coulomb interactions. This equation is computationally unviable for all but the 
smallest molecules owing to the exponential barrier of complexity with the number of electrons 
m- Therefore, approximated models have been developed to study such systems. 


3.1. Kohn-Sham formulation 

In the early 60’s Kohn and Hohenberg proved that the ground-state energy of a quantum many- 
body system is a unique functional of the one-body electron density p(f) [2]. Starting from 
this basic result, shortly later Kohn and Sham, proceeded to derive a set of effective one-body 
time-independent Schrodinger equations (SE) of the form [3]: 


Hxsfij = Ej(f>j , 


( 2 ) 


where is the wave function of the j-th orbital with energy Ej and the KS Hamiltonian 

reads as follows: 


Hrs 


h 2 2 
_ 2m Vj , 


P(f 


\r — r 


-d\' + y. z i 

i 


p(r') 

I Ri-r* 


-d \' + V x , 


[p] 


( 3 ) 


where Rj and Zj are the coordinates and charges of the nuclei, respectively. In the above, V xc 
defines the exchange-correlation energy, the term collecting all of the unknown effects of many- 
body interactions and Pauli exclusion principle. Note that, by effect of the Kohn-Hohenberg 
theorem, this term is guaranteed to depend only on the total electron density 

p(r) = J2 l<M f )| 2 > ( 4 ) 

j 

which represents an enormous simplification of the quantum many-body problem. If such 
dependence could be sorted out, the Kohn-Sham equations (KSE) would be exact. 

From the mathematical point of view, the KSE’s represent a set of single-particle time- 
independent Schrodinger equations, coupled only through the total density p(r). Furthermore, 
the orbitals must obey the orthogonality condition 

J cfij (r) fik if) dr = djh , 


( 5 ) 




which is a set of N e (N e + l)/2 global constraints. To be noted that the above orbitals are a 
mere computational device to perform variational minimisation of the many-body ground-state 
energy, with no specific physical meaning. 


4. The kinetic approach to Kohn-Sham theory 

The main observation is that, upon performing a Wick rotation, t t' = —it, the time- 
dependent Schrodinger equation is basically an one-particle diffusion-reaction equation, where 
the reaction term contains all the details of the potential energy, i.e., 
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where ^ is the wave function, and V is the potential energy. Since reaction-diffusion equations 
are known to emerge as a macroscopic limit of an underlying Boltzmann equation, it is natural 
to treat the wave function as a “fluid” and propose a (lattice) kinetic representation for the wave 
function, such that 

^Zfp{r-,t') =ip(r;t') , (7) 

p 

In the above, / p (r), with p = 0, is the probability of finding a “particle” with velocity 
v p at position f and time t, where b is the total number of discrete velocity vectors. For the 
moment, “particle” means simply a computational quasi-particle, with no physical implications 
for the existence of an underlying kinetic theory behind the Schrodinger representation. Note 
that assuming an arbitrary initial condition for the wave function if, one can expand it using a 
basis of orthogonal orbitals, 


iP(r,t') = ^2 a k (pk(r) exp(-E k t'/h) 
k 


( 8 ) 


such that after some time, it will be ijj(f,t') ~ no0o(O exp(— E^t'/fi), which upon normalisation 
leads to the ground state solution of the time-independent Schrodinger equation, 
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The distribution function f p evolves according to the lattice Boltzmann equation (time step 
made unity for simplicity), 


f P {r + v p ;t + l) = f p (f;t)-u}(f p -fp) + V p , (10) 

where w is a typical relaxation frequency which controls the quantum diffusivity D = Ti/m , 
fp is a suitable local equilibrium ensuring mass conservation, and Ip is a suitable source term 
representing the particle loss/gain associated with the potential energy scattering. In lattice 
units Ax = At = 1, one has: 

D = cl( 1/a; - 1/2) , (11) 

where c s is the lattice sound speed. Thus, by changing ui, one has a handle on the effective mass 
of the electrons. The explicit expression for the equilibrium distribution is as follows: 


4 e = Hpip , (12) 

where = (1 + (D — c 2 )q p ), with w p being suitable lattice weights normalised to unity, and 
q p = ( v 2 — 3c^)/2c^ the second order lattice Hermite polynomial. The scattering source V p reads 
as follows: 

V p = x.pji’ , 


( 13 ) 


where Xp = Wp[ 1 — (u^ — 3c^)/2c^)] ensures that the first and second order moments are zero, 
namely J2 P fpVp = 0 and Yh P fp^pVp = 0, so that momentum and energy of the wave function 
are unaffected. 

In order to solve the KS system, we need a distribution function for each jth-orbital, i.e. 
fp -4 fjp and consequently i/j — > i/jj, and the potential energy V contains the electron-electron 
and ion-electron interactions, plus the exchange-correlation potential. Thus, the corresponding 
electronic Lattice Boltzmann (LB) equation reads as follows: 


f jp (r + v p ; t + 1) = f jp (r-1) - u{f jp - f- p ) + V jp + W jp , (14) 

where, Wj p is the source term in charge of securing the orthogonality constraints (j5|). The 
explicit expressions for the equilibrium distributions and scattering potentials are as follows: 


V 

fjp = Zp'l’j, and v iP = X-Pj-'fj 


and the orthogonalisation potential takes the form 
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where A j k = fi’jl'^k )/is the cosine of the angle between orbitals j and k in Hilbert space. 
Clearly, once two orbitals are mutually orthogonal, A jk = 0, so that they no longer contribute 
to the potential Wj p . 

In order to understand the role of the orthogonalisation potential we will introduce the 
following derivation. Let us first assume that all orbitals are initialised with the same wave 
function, ipj(r) = ^(r), and we do not consider any orthogonalisation potential in equation (114[) . 
i.e. Wj p = 0. According to our previous explanation, the wave functions can be written as 
expansions in the orthogonal basis ^.(r), 


ipj(r, t') = Y] fjp = Y a kj<l>k(r) exp(-E k t'/h) . (17) 

P k 

where a k j = {^kl'f’j) are projection coefficients. After some time, the only leading term in the 
sum will be the ground state orbital, </>o, such that ^j(f,t') ~ aoj(/)o(r) exp(— Eot'/h). At this 
stage, all orbitals in our lattice kinetic approach will reach the same ground state. Let us now 
subtract from the wave function ipi, the wave function i/jq, such that 

ipi(r,t') = Yj a kif>k(f) exp(—E k l//ti) - AioV’o {r,t') . (18) 

k 

This equation can be rewritten explicitly as, 

= Y ~ (<fo#o)) 4>k(r) e*p(-E k t'/h) . (19) 

The evolution equation for the orbital ipo implies that it will converge eventually to the ground 
state (j>Q. Therefore, the first term in the sum, for k = 0, vanishes and the remaining contributions 
to the sum are 

= Y a ki4>k(r) exp(-E k i'/h) 

fc >0 


( 20 ) 



Thus, the orbital 'ipi will eventually converge to the first excited state <f>\. Repeating this 
procedure for each wave function, one can write 
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or 

-0 j(r,t ') = Y a kj^k{f) exp (~E k t'/h) - Y , (22) 

k Kj 

and show that each wave function will eventually converge, upon normalisation, to the respective 
KS orbital, i.e. ipj —> cf)j. In order to include this procedure into the evolution of the wave 
functions, we need to force the jth wave function to converge to 

i>j if, £') = ipj (r, t')~Y (r, t ') = Y fop ~ T (r, t') . (23) 

Kj v Kj 


This can be done by replacing in the equilibrium distribution ijjj by if)*-, obtaining the following 
lattice kinetic equation, 


fj p {r + v p ; t + 1) = f jp (f] t) - u(f jp - £pipj) + V jp . (24) 

By reorganising terms, this equation leads to 

fjp( r + v p'i t T 1) = fjp{ r 'i t) — &(fjp — £ p ipj) + Vjp £ p uj 'y ] A t ) , (25) 

Kj 

which is precisely our previous definition of the lattice kinetic equation with the orthogonalisation 
potential. Note that after some time, we have ipj —> 0| —> This process is a kind of time- 
dependent Gram-Schmidt procedure. Indeed, one can also evolve first the wave function 0 O , 
and once the ground state is found, evolve 0i, and so on. However, this will lead to longer 
computational times to achieve the total electronic ground state. The scaling properties and 
convergency speed of the orthogonalisation procedure described here will be a subject of future 
research. 


4-.1. Concurrent Dynamics 

The above electronic LB only applies to ground-state calculations, where the electronic 
distribution is systematically adjusted to the Born-Oppenheimer surface defined by the actual 
positions of the ions, Ri(t), these latter being obtained by a standard Molecular Dynamics 
integration. The total energy is thus evolving in (imaginary) time until the total ground state 
value is attained. This is the standard Born-Oppenheimer (BO) scenario. 

The KS-LB can be extended to the concurrent-dynamic (CD) scenario, whereby the electrons 
and the ions are moved simultaneously, according to the celebrated Car-Parrinello picture p3|- In 
this scenario, energy owes to be conserved, and consequently the irreversible relaxation operator 
in Eq. (|14l) needs to be augmented with a potential, U, such that density and momentum are 
conserved, while energy is supplied to (removed from) the system so as to keep it constant. 

The corresponding LB equation takes then the form: 

fjp(r + v p h ; t + h) = fj p (f-v p h ; t - h) + 2 h(U jp + iV jp + iW jp ) , (26) 

where h is the time-step. The specific expression of Uj p is as follows: 


Ujp — KWp Qp Ej (f)j , 


( 27 ) 





Figure 1. Convergence of the total energy for the atoms H and Be as a function of the number 
of iterations performed by our model. The dashed-line corresponds to the final value, E = 0.500 
a.u. (1 a.u. ~ 27.2 eV) and E = 14.45 a.u., for the hydrogen and beryllium atoms, respectively. 
For Beryllium, we found an error of 1% with respect to the expected value in Refs. [15 , om, 
while for hydrogen the error is of the order of 0.1%. We have taken tk = St. 


Atom 

v x 

Exp. V x 

Vc 

Exp. V c 

Time 

Iterations 

H 

-0.310 

-0.310 

0 

0 

9 sec 

158 

He 

-1.025 

-1.025 

-0.044 

-0.044 

1 min 

906 

Be 

-2.658 

-2.658 

-0.095 

-0.095 

3 min 

965 

B 

-3.728 

-3.728 

-0.128 

-0.128 

29 min 

2704 

C 

-5.032 

-5.032 

-0.161 

-0.161 

33 min 

2305 


Table 1. Exchange-Correlation energies for H, He, Be, B, and C in atomic units (a.u.). 
Computational time and the number of iterations performed by the model to reach the ground 
state are also shown. The expected values of exchange, V x , and correlation, V c , energies are 
taken from Refs. Hi- The simulations have run on an Intel Core i7 at 2.3 GHz. We have taken 
tk = St. 


k being an adjustable parameter in charge of keeping the energy at Ej. To be noted that the 
algorithm is now a two-step time-centered scheme, in compliance with the requirement of energy 
conservation. This is the analogue, in kinetic language, of Car-Parrinello’s inertial term in the 
electron equation, i.e. a second order time derivative to comply with the overarching electron-ion 
Lagrangian. The fact that the present LB scheme operates successfully also in the CD regime 
hints at the existence of an overarching kinetic Lagrangian. At the time of this writing, however, 
such Lagrangian remains to be found. 

5. Results 

In order to show the validity of our approach, we present the calculations of the exchange and 
correlation energies for different atoms, H, He, Be, B, and C. For H, Be, B, and C, one can 
find more details about the calculations and parameter values in Ref. m- For those cases, 
we have used lattice sizes of 32 3 , 40 3 , 56 3 , and 58 3 , respectively. For He, which has not been 
reported before using our model, we have taken a system size of 36 3 , Bohr radius a q = 6.5, 
Tl/m = 1, and r = 1 (all values are in numerical units). We have also used fourth order in 
the Hermite expansion for the equilibrium distribution and source term (for details about the 





















Figure 2. Water molecule, H 2 O. The blue and red isosurfaces denote low and high electronic 
density, respectively. Using our model, we have obtained for the angles between bonds, 104.4 
degrees, and a O-H bond distance of 0.95 A. 


Hermite expansion, see Ref. [12]). The measured values for the exchange and correlation energies 
are given in Table [U together with the computational time and number of iterations. Note that 
the reported number of iterations might exceed the one of existing iterative methods. However, 
as we shall show shortly, each of our iterations takes typically much less computational time, thus 
making our method competitive. In Fig. [TJ we see that the convergence of the total energy, for 
H and Be, to the ground state is very fast in the beginning and slows down when it approaches 
the exact solution. The differences between the obtained and expected exchange and correlation 
energies are of the order of 1%. The simulations are stopped when the energy of the orbitals 
presents changes of less than 10~ 5 % between two subsequent steps. 

As an additional application, we build the water molecule from scratch, by using the bare 
Coulomb potential. For this simulation, we use a lattice size of 74 3 , %/m = 1, ao = 10.5 cells, and 
M = 5, and place the oxygen atom at the center of the lattice. The hydrogen atoms are located 
randomly in space, and we let the system evolve to the ground configuration. After 4 hours 
(10447 iterations), we achieved the configuration shown in Fig. [2j The angle, 104.4 degrees, and 
the bond distance of 0.95 A are in excellent agreement with expected and experimental values 

mm- 

Finally, we compare BO and CD versions of our kinetic scheme against each other. Here, we 
assume that the BO dynamics provides the correct results, and will be used to validate the CD. 
To this purpose, we excite the H 2 molecule and let it vibrate in its first mode. We have used 
a lattice size of 24 3 . Further details about the values of the model parameters can be found 
in Ref. [12]. For the BO dynamics we use the equation of motion for the ions and, at each 
time step, we calculate the ground state of the electronic density using Eq. (| 14[) . On the other 
hand, for CD we evolve simultaneously the ions and the electronic orbitals, using the discrete 
version of Eq. (126|) . Since in CD there is no relaxation process to the ground state, we must use 
a smaller time step such that electrons can accommodate according to the new position of the 
ions at each time step. The time step of the BO molecular dynamics is set to At bo = 0-09 fs, 
while for CD, AtcD = 0.0014 fs. Although the time step in CD is smaller by nearly two orders 



t(fs) 


Figure 3. Vibrational mode of the diatomic molecule H 2 . Here BO and CD denote Born- 
Oppenheimer and Concurrent Dynamics, respectively. The inset shows that up to 500 fs, the 
amplitude of the oscillations does not show any appreciable decay in time. Here, ao is the Bohr 
radius, dn-H the distance between H atoms and d g ^_ H the ground state distance. 


of magnitude, it runs 17 times faster (taking 17 sec to simulate one fs) than BO dynamics. In 
fig. H we report the first vibrational mode of the hydrogen molecule with BO and CD after 40 
fs, from which one can appreciate that the trajectories of CD are pretty close to the respective 
BO ones. In the inset of this figure, we can also see that after 500 fs, the amplitude of the 
oscillations remains constant, without any appreciable loss of energy. The oscillation frequency 
of the vibrational mode obtained by the simulation is 4206 cm” 1 , which departs by 1% from the 
experimental value of 4161 cm” 1 |19l . 

Interestingly, the LB simulations do resolve the full Coulomb interactions, without making 
any use of soft pseudo-potentials pj2]. This is possibly due to the high-order of spatial accuracy 
achieved by the source terms. The amount of computational memory that our model uses can 
be roughly estimated as Mt = 8(36 + 2No)V in bytes, where No is the number of orbitals and 
V the lattice size. For instance, for the water molecule the total memory used is about 150 
MB, and by using a standard computer (or a single node of a cluster) with 8 GB of RAM, one 
could study about 8 of these molecules. This small system size implies that parallel computing 
is a crucial step in order to study larger system sizes, and the above relation for the memory 
requirements provides a reasonable estimate of the maximum system size that one can simulate 
at a given level of computing power. 

Our simulations have been performed using open boundary conditions. However, the use 
of periodic boundary conditions is also possible at the price of introducing corrections due to 
Ewald summation for the electric potential. Indeed, using periodic boundary conditions can lead 
to more stable simulations, since the distribution functions do not need to be extrapolated at 
the boundaries, and therefore, their dynamics or relaxation to ground state follows the correct 
equations holding in the bulk. 

6. Prospective developments 

The kinetic representation discussed in this paper lends itself to a number of prospective 
developments. First, one should note that such representation is not restricted to the time- 
independent Kohn-Shan formulation. Indeed, in the recent past Runge and Gross have extended 
Kohn-Hohenberg’s theorem to the case of time-dependent potentials |2Q]. In this case, the KS 
formalism still applies, although in real-time form, i.e: 


ikd t (j)j(r;t ) = H KS (t)<f>j(r;t ) 


( 28 ) 





























where the KS Hamiltonian now contains a time-dependent potential, Hks = K + V(r;t). 

The above equations can be integrated using consolidated quantum LB schemes m, 
exporting the very same upgrades described in this paper. Another direction, possibly even 
closer in spirit to kinetic theory, is the so-called density-current formulation of TDDFT, in which 
the hamiltonian depends not only on the electron density, but also on the electron current: 

= Y^^i\M ^ t )\ 2 ( 29 ) 
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where vj = is the velocity operator. 

The density and current functionals evolve according to the continuity equation: 

d t p + V • J = 0 (30) 

which provides p(r\t ) once the current functional J = J(p(f-,t )) is prescribed. Of course, 
the latter is not known exactly, but sensible approximations have been developed in the 
specialised literature, e.g the Vignale-Kohn functional m, which could be implemented in 
the corresponding electronic LB exactly along the lines described in this paper. 

Finally, one may also envisage a more radical strategy, whereby the energy functional would 
be expressed in terms of the Boltzmann distribution / itself rather than by its lower-order kinetic 
moments, such as density, current, and momentum flux. The prospective Boltzmann functional 
would symbolically read as 

E = j £[f]dfdv (31) 

where the double integration extends over the six-dimensional phase space. 

What would one gain out of such functional? Possibly, the non-perturbative incorporation of 
inhomogeneity effects at all orders. Indeed, the Boltzmann distribution splits naturally into an 
equilibrium and non-equilibrium component, / = f e + f ne , the former depending parametrically 
only on invariants, density and momentum, while the latter, by definition, collects the effects of 
space-time inhomogeneity at all orders. Of course, such functional would make the difference 
only at the level where further moments on top of current-density, are required. 

7. Conclusions 

Summarising, we have presented the main ideas behind the kinetic formulation of electronic 
density functional theory, along with the ensuing implementation within the Lattice Boltzmann 
formalism. We have described in detail the effects of the orthogonalisation potential and shown 
how it contributes to find the system of orthogonal KS orbitals. In order to show the validity 
of the lattice Boltzmann approach for KS, we have shown numerical data for the computation 
of the exchange and correlation energies for simple atoms, finding excellent agreement with 
the previous literature. The results of numerical simulations for the first vibrational mode of 
the hydrogen molecule are also presented. The prospects, including future extensions to time- 
dependent DFT, look exciting. 
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